* program returning potential (0) or discharge (1) for
* steady (0) or transient step (1) quadrilateral areasink
* to be run from Python program with
*    import popen2 // pipe open 2 (1 for input + 1 for output)
*    fin,fout=popen2.popen2('SAP0')
*    fout.write('pot0dis1 st0tns1\n x y x1 y1 x2 y2 x3 y3 x4 y4\n')
*    a=fin.readline()
*    a=eval(a)
* where a is the value of the influence function, or above plus
*    b=fin.readline()
*    b=eval(b)
* where XY(a,b) is the value of the influence function
      program SApq01
      implicit none
      logical          testvar
      integer          ipot0dis1, ist0tns1
      double precision x, y, dt,
     &                 x1, y1, x2, y2, x3, y3, x4, y4, dp,
     &                 RFSAP0, RFSAP1
      complex*16       cz, cz1, cz2, cz3, cz4, cq,
     &                 CFSAQ0, CFSAQ1
      logical          COMMITQQ

      testvar=COMMITQQ(6)
      do while(.true.)
         read(5,*,err=9000) ipot0dis1, ist0tns1
         if (ipot0dis1.ne.0 .and. ipot0dis1.ne.1) then
*           neither potential nor discharge: QUIT
            goto 9000
         end if
         if (ist0tns1.eq.0) then
*           steady
            read(5,*,err=100,end=9000) x, y,
     &         x1, y1, x2, y2, x3, y3, x4, y4
         else if (ist0tns1.eq.1) then
*           transient step
            read(5,*,err=100,end=9000) x, y, dt,
     &         x1, y1, x2, y2, x3, y3, x4, y4, dp
         else
*           neither steady nor transient step: ERROR
            goto 100
         end if

         cz  = cmplx(x ,y )
         cz1 = cmplx(x1,y1)
         cz2 = cmplx(x2,y2)
         cz3 = cmplx(x3,y3)
         cz4 = cmplx(x4,y4)

         if (ipot0dis1.eq.0) then

*           potential
            if (ist0tns1.eq.0) then
*              steady
               write(6,*) RFSAP0(cz,  cz1, cz2, cz3, cz4)
            else if (ist0tns1.eq.1) then
*              transient step
               write(6,*) RFSAP1(cz, dt,  cz1, cz2, cz3, cz4,  dp)
            end if

         else if (ipot0dis1.eq.1) then

*           discharge
            if (ist0tns1.eq.0) then
*              steady
               cq  = CFSAQ0(cz,  cz1, cz2, cz3, cz4)
            else if (ist0tns1.eq.1) then
*              transient step
               cq  = CFSAQ1(cz, dt,  cz1, cz2, cz3, cz4,  dp)

            end if
            write(6,*)  real(cq)
            write(6,*) aimag(cq)

         end if

         goto 200

100      write(6,'(''1.e32'')')
         if (ipot0dis1.eq.1) write(6,'(''1.e32'')')

200      continue
         call flush(6)
      end do

9000  end program
